|
|
(******************************************************************************)
(**) ОТДЕЛ Матр;
(******************************************************************************
* НАЗНАЧЕНИЕ: определение матриц и основные вычисления над ними
*
* Источники, ссылки, библиография
Debord J., Библиотека математических подпрограмм
tpmath.zip
Goffe B., Программа SimAnn.FOR (Simulated Annealing in Fortran)
http://www.netlib.org/opt/simann.f
EISPACK: Библиотека Фортран функций для вычисления собственных значений и векторов
http://www.netlib.org/eispack
Marsaglia G., Тесты для генераторов случайных чисел
http://stat.fsu.edu/~geo/diehard.html
Moshier S., Библиотека математических подпрограмм
http://www.moshier.net
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.
Численные рецепты на Паскале
http://garbo.uwasa.fi/pc/turbopas.html
Пакет численных методов для Турбо Паскаля
Borland International, 1986
*
******************************************************************************)
ИСПОЛЬЗУЕТ Матем,Вект;
ВИД
Вещ = Матем.Вещ;
(******************************************************************************
* Код ошибки после выполнения функций
******************************************************************************)
ПОСТ
ВЫРОЖДЕНА- =-1; (* вырожденная матрица *)
НЕТ_СХОДИМОСТИ-=-2; (* нет сходимости *)
НЕ_П_О- =-3; (* матрица не является положительно определённой *)
ВИД
Вид-=РЯД ИЗ Вект.Вид;
Доступ-=ДОСТУП К Вид;
Перестановки-=ДОСТУП К РЯД ИЗ ЦЕЛ;
(******************************************************************************
* Переписывание матриц и векторов
******************************************************************************)
(* заготовки ЗАДАЧ *)
ЗАДАЧА^ ОбменСтрок-(A+:Вид; i,k:ЦЕЛ);
ЗАДАЧА^ ОбменСтолбцов-(A+:Вид; j,k:ЦЕЛ);
ЗАДАЧА^ СписатьВектор-(из-,в+:Вект.Вид);
ЗАДАЧА^ Списать-(из-,в+:Вид);
ЗАДАЧА^ СписатьСтолбецИзВектора-(из-:Вект.Вид; в+:Вид; стлб:ЦЕЛ);
ЗАДАЧА^ СписатьВекторИзСтолбеца-(из-:Вид; в+:Вект.Вид; стлб:ЦЕЛ);
(******************************************************************************)
ЗАДАЧА^ МетодГаусса-(A-:Вид; b-:Вект.Вид; Ao+:Вид; x+:Вект.Вид):ЦЕЛ;
(* Цель: решение системы линейных алгебраических уравнений методом
* исключения Гаусса
******************************************************************************
* До: < A > - матрица системы
* < b > - вектор свободных членов
* После: < Ao > - обратная матрица
* < x > - вектор решения
* Ответ: 0 или ВЫРОЖДЕНА *)
(******************************************************************************)
ЗАДАЧА^ Обратить-(A-,Ao+:Вид):ЦЕЛ;
(* Цель: обращение квадратной матрицы методом Гаусса
******************************************************************************
* До: < A > - исходная матрица
* После: < Ao > - обратная матрица
* Ответ: 0 или ВЫРОЖДЕНА *)
(******************************************************************************)
ЗАДАЧА^ Определитель-(A+:Вид):Вещ;
(* Цель: вычисление определителя квадратной матрицы
******************************************************************************
* До: < A > - исходная матрица
* После: < A > - разрушенная матрица
* Ответ: значение определителя *)
(******************************************************************************)
ЗАДАЧА^ РазложитьНаLL-(A-,L+:Вид):ЦЕЛ;
(* Цель: Разложение Холецкого (A=L*L') симметрической положительно
* определенной матрицы A на две треугольные матрицы L и L', где L - это
* нижняя треугольная матрица, а L' - матрица, транспонированная к L.
******************************************************************************
* Прим: Эту задачу можно использовать для испытания на положительно
* определенность.
* До: < A > - исходная матрица A
* После: < L > - матрица L
* Ответ: 0 или НЕ_П_О *)
(******************************************************************************)
ЗАДАЧА^ РазложитьНаLU-(A+:Вид; строки+:Перестановки):ЦЕЛ;
(* Цель: LU-разложение (A=L*U) квадратной матрицы A на два сомножителя L и U,
* где L - это нижняя треугольная матрица с единичной главной диагональю,
* а U - верхняя треугольная матрица. Эта задача используется
* совместно с РешитьИзLU для решения систем уравнений.
******************************************************************************
* До: < A > - исходная матрица
* После: < A > - содержит элементы и L и U
* <строки> - вновь созданные записи о перестановках строк
* Ответ: 0 или ВЫРОЖДЕНА *)
(******************************************************************************)
ЗАДАЧА^ РешитьИзLU-(A-:Вид; b-:Вект.Вид; строки:Перестановки; x+:Вект.Вид);
(* Цель: решение системы уравнений после преобразования исходной матрицы
* с помощью РазложитьНаLU
* До: < A > - матрицы L и U после РазложитьНаLU
* < b > - вектор свободных членов
* <строки> - записи о перестановках строк из РазложитьНаLU
* После: < x > - вектор решения *)
(******************************************************************************)
ЗАДАЧА^ РазложитьНаSV-(A+:Вид; s+:Вект.Вид; V+:Вид):ЦЕЛ;
(* Цель: SV разложение прямоугольной матрицы A (N x M, с N >= M) на три
* сомножителя (A=U*S*V'), где U (N x M) состоит из столбцов - собственных
* векторов, S (M x M) - диагональная матрица с элементами >= 0 (собственные
* значения) и V (M x M) - ортогональная матрица. Эта задача используется
* совместно с РешитьИзSV для решения систем уравнений.
******************************************************************************
* До: < A > - исходная матрица
* После: < A > - содержит элементы U
* < s > - вектор собственных значений
* < V > - ортогональная матрица
* Ответ: 0 или НЕТ_СХОДИМОСТИ *)
(******************************************************************************)
ЗАДАЧА^ ОбнулитьS-(s+:Вект.Вид; откл:Вещ);
(* Цель: обнуляет наименьшие собственные значения по заданной границе
******************************************************************************
* До: < s > - вектор собственных значений S
* <откл> - относительное отклонение (граница будет равна откл*макс(S) )
* После: < s > - измененный вектор собственных значений *)
(******************************************************************************)
ЗАДАЧА^ РешитьИзSV-(U-:Вид; s-:Вект.Вид; V-:Вид; b-,x+:Вект.Вид);
(* Цель: решение системы уравнений после преобразования исходной матрицы
* с помощью РазложитьНаSV и обнуления наименьших собственных значений
* с помощью ОбнулитьS
******************************************************************************
* До: < U > - матрица U после РазложитьНаSV
* < s > - вектор S после РазложитьНаSV
* < V > - матрица V после РазложитьНаSV
* < b > - вектор свободных членов
* После: < x > - вектор решения (V*Diag(1/S(i))*U'*B, для S(i) # 0) *)
(******************************************************************************)
ЗАДАЧА^ ВосстановитьИзSV-(U-:Вид; s-:Вект.Вид; V-,A+:Вид);
(* Цель: восстанавливает матрицу A перемножая U, S и V' после
* обнуления наименьших собственных значений задачей ОбнулитьS
******************************************************************************
* До: < U > - матрица U после РазложитьНаSV
* < s > - вектор S после РазложитьНаSV
* < V > - матрица V после РазложитьНаSV
* После: < A > - восстановленная матрица *)
(******************************************************************************)
ЗАДАЧА^ РазложитьНаQR-(A+,R+:Вид):ЦЕЛ;
(* Цель: QR разложение прямоугольной матрицы (N x M, с N >= M) на два
* сомножителя (A=Q*R), где Q (N x M) - ортогональная матрица, а
* R (M x M) - верхняя треугольная матрица. Эта задача используется
* совместно с РешитьИзQR для решения систем уравнений.
******************************************************************************
* До: < A > - исходная матрица
* После: < A > - содержит элементы Q
* < R > - верхняя треугольная матрица
* Ответ: 0 или ВЫРОЖДЕНА *)
(******************************************************************************)
ЗАДАЧА^ РешитьИзQR-(q-,R-:Вид; b-,x+:Вект.Вид);
(* Цель: решение системы уравнений после преобразования исходной матрицы
* с помощью РазложитьНаQR
******************************************************************************
* До: < q > - матрица Q после РазложитьНаQR
* < R > - матрица R после РазложитьНаQR
* < b > - вектор свободных членов
* После: - вектор решения *)
(* полные реализации ЗАДАЧ *)
(******************************************************************************)
ЗАДАЧА ОбменСтрок-(A+:Вид; i,k:ЦЕЛ);
ПЕР
j:ЦЕЛ;
УКАЗ
ОТ j:=0 ДО РАЗМЕР(A,1)-1 ВЫП
Матем.обмен(A[i,j],A[k,j])
КОН
КОН ОбменСтрок;
(******************************************************************************)
ЗАДАЧА ОбменСтолбцов-(A+:Вид; j,k:ЦЕЛ);
ПЕР
i:ЦЕЛ;
УКАЗ
ОТ i:=0 ДО РАЗМЕР(A,0)-1 ВЫП
Матем.обмен(A[i,j],A[i,k])
КОН
КОН ОбменСтолбцов;
(******************************************************************************)
ЗАДАЧА СписатьВектор-(из-,в+:Вект.Вид);
ПЕР
i:ЦЕЛ;
УКАЗ
ОТ i:=0 ДО РАЗМЕР(из)-1 ВЫП
в[i]:=из[i]
КОН
КОН СписатьВектор;
(******************************************************************************)
ЗАДАЧА Списать-(из-,в+:Вид);
ПЕР
i,j:ЦЕЛ;
УКАЗ
ОТ i:=0 ДО Матем.минЦ(РАЗМЕР(из,0),РАЗМЕР(в,0))-1 ВЫП
ОТ j:=0 ДО Матем.минЦ(РАЗМЕР(из,1),РАЗМЕР(в,1))-1 ВЫП
в[i,j]:=из[i,j]
КОН
КОН
КОН Списать;
(******************************************************************************)
ЗАДАЧА СписатьСтолбецИзВектора-(из-:Вект.Вид; в+:Вид; стлб:ЦЕЛ);
ПЕР
i:ЦЕЛ;
УКАЗ
ОТ i:=0 ДО Матем.минЦ(РАЗМЕР(из),РАЗМЕР(в,0))-1 ВЫП
в[i,стлб]:=из[i]
КОН
КОН СписатьСтолбецИзВектора;
(******************************************************************************)
ЗАДАЧА СписатьВекторИзСтолбеца-(из-:Вид; в+:Вект.Вид; стлб:ЦЕЛ);
ПЕР
i:ЦЕЛ;
УКАЗ
ОТ i:=0 ДО Матем.минЦ(РАЗМЕР(из),РАЗМЕР(в,0))-1 ВЫП
в[i]:=из[i,стлб]
КОН
КОН СписатьВекторИзСтолбеца;
(******************************************************************************)
ЗАДАЧА МетодГаусса-(A-:Вид; b-:Вект.Вид; Ao+:Вид; x+:Вект.Вид):ЦЕЛ;
(* Цель: решение системы линейных алгебраических уравнений методом
* исключения Гаусса
******************************************************************************
* До: < A > - матрица системы
* < b > - вектор свободных членов
* После: < Ao > - обратная матрица
* < x > - вектор решения
* Ответ: 0 или ВЫРОЖДЕНА *)
ПЕР
i,j,k,посл:ЦЕЛ;
вед,t:Вещ;
ведСтр,ведСтлб:ДОСТУП К РЯД ИЗ ЦЕЛ; (* строки и столбцы ведущего элемента *)
УКАЗ
посл:=РАЗМЕР(A)-1;
СОЗДАТЬ(ведСтр,посл+1);
СОЗДАТЬ(ведСтлб,посл+1);
Списать(A,Ao);
СписатьВектор(b,x);
ОТ k:=0 ДО посл ВЫП
(* поиск наибольшего ведущего элемента в подматрице Ao[k..посл,k..посл] *)
вед:=Ao[k,k];
ведСтр[k]:=k;
ведСтлб[k]:=k;
ОТ i:=k ДО посл ВЫП
ОТ j:=k ДО посл ВЫП
ЕСЛИ МОДУЛЬ(Ao[i,j]) > МОДУЛЬ(вед) ТО
вед:=Ao[i,j];
ведСтр[k]:=i;
ведСтлб[k]:=j
КОН
КОН
КОН;
(* ведущий элемент слабоват => псевдо-вырожденная матрица *)
ЕСЛИ МОДУЛЬ(вед) < Матем.МАШЕПС ТО
ВОЗВРАТ ВЫРОЖДЕНА
КОН;
(* обмен текущей строки (k) со строкой ведущего элемента *)
ЕСЛИ ведСтр[k] # k ТО
ОбменСтрок(Ao,ведСтр[k],k);
Матем.обмен(x[ведСтр[k]],x[k])
КОН;
(* обмен текущего столбца (k) со столбцом ведущего элемента *)
ЕСЛИ ведСтлб[k] # k ТО
ОбменСтолбцов(Ao,ведСтлб[k],k)
КОН;
(* преобразование строки ведущего элемента *)
Ao[k,k]:=1;
ОТ j:=0 ДО посл ВЫП
Ao[k,j]:=Ao[k,j]/вед
КОН;
x[k]:=x[k]/вед;
(* преобразование оставшихся строк *)
ОТ i:=0 ДО посл ВЫП
ЕСЛИ i # k ТО
t:=Ao[i,k];
Ao[i,k]:=0;
ОТ j:=0 ДО посл ВЫП
Ao[i,j]:=Ao[i,j]-t*Ao[k,j]
КОН;
x[i]:=x[i]-t*x[k]
КОН
КОН
КОН;
(* перестановка обратной матрицы *)
ОТ i:=посл ДО 0 ПО -1 ВЫП
ЕСЛИ ведСтлб[i] # i ТО
ОбменСтрок(Ao,ведСтлб[i],i);
Матем.обмен(x[ведСтлб[i]],x[i])
КОН
КОН;
ОТ j:=посл ДО 0 ПО -1 ВЫП
ЕСЛИ ведСтр[j] # j ТО
ОбменСтолбцов(Ao,ведСтр[j],j)
КОН
КОН;
ВОЗВРАТ 0
КОН МетодГаусса;
(******************************************************************************)
ЗАДАЧА Обратить-(A-:Вид; Ao+:Вид):ЦЕЛ;
(* Цель: обращение квадратной матрицы методом Гаусса
******************************************************************************
* До: < A > - исходная матрица
* После: < Ao > - обратная матрица
* Ответ: 0 или ВЫРОЖДЕНА *)
ПЕР
i,j,k,посл:ЦЕЛ;
вед,t:Вещ;
ведСтр,ведСтлб:ДОСТУП К РЯД ИЗ ЦЕЛ; (* строки и столбцы ведущего элемента *)
УКАЗ
посл:=РАЗМЕР(A)-1;
СОЗДАТЬ(ведСтр,посл+1);
СОЗДАТЬ(ведСтлб,посл+1);
Списать(A,Ao);
ОТ k:=0 ДО посл ВЫП
(* поиск наибольшего ведущего элемента в подматрице Ao[k..посл,k..посл] *)
вед:=Ao[k,k];
ведСтр[k]:=k;
ведСтлб[k]:=k;
ОТ i:=k ДО посл ВЫП
ОТ j:=k ДО посл ВЫП
ЕСЛИ МОДУЛЬ(Ao[i,j]) > МОДУЛЬ(вед) ТО
вед:=Ao[i,j];
ведСтр[k]:=i;
ведСтлб[k]:=j
КОН
КОН
КОН;
(* ведущий элемент слабоват => псевдо-вырожденная матрица *)
ЕСЛИ МОДУЛЬ(вед) < Матем.МАШЕПС ТО
ВОЗВРАТ ВЫРОЖДЕНА
КОН;
(* обмен текущей строки (k) со строкой ведущего элемента *)
ЕСЛИ ведСтр[k] # k ТО
ОбменСтрок(Ao,ведСтр[k],k)
КОН;
(* обмен текущего столбца (k) со столбцом ведущего элемента *)
ЕСЛИ ведСтлб[k] # k ТО
ОбменСтолбцов(Ao,ведСтлб[k],k)
КОН;
(* преобразование строки ведущего элемента *)
Ao[k,k]:=1;
ОТ j:=0 ДО посл ВЫП
Ao[k,j]:=Ao[k,j]/вед
КОН;
(* преобразование оставшихся строк *)
ОТ i:=0 ДО посл ВЫП
ЕСЛИ i # k ТО
t:=Ao[i,k];
Ao[i,k]:=0;
ОТ j:=0 ДО посл ВЫП
Ao[i,j]:=Ao[i,j]-t*Ao[k,j]
КОН
КОН
КОН
КОН;
(* перестановка обратной матрицы *)
ОТ i:=посл ДО 0 ПО -1 ВЫП
ЕСЛИ ведСтлб[i] # i ТО
ОбменСтрок(Ao,ведСтлб[i],i)
КОН
КОН;
ОТ j:=посл ДО 0 ПО -1 ВЫП
ЕСЛИ ведСтр[j] # j ТО
ОбменСтолбцов(Ao,ведСтр[j],j)
КОН
КОН;
ВОЗВРАТ 0
КОН Обратить;
(******************************************************************************)
ЗАДАЧА Определитель-(A+:Вид):Вещ;
(* Цель: вычисление определителя квадратной матрицы
******************************************************************************
* До: < A > - исходная матрица
* После: < A > - разрушенная матрица
* Ответ: значение определителя *)
ПЕР
d,t:Вещ; (* частичный определитель и множитель *)
i,j,k,посл:ЦЕЛ;
опред0:КЛЮЧ; (* если определитель = 0 *)
УКАЗ
посл:=РАЗМЕР(A)-1;
опред0:=ОТКЛ;
d:=1;
(* строим верхнюю треугольную матрицу *)
k:=0;
ПОКА НЕ опред0 И (k < посл) ВЫП
(* при 0-м диагональном элементе переставляем строки *)
ЕСЛИ МОДУЛЬ(A[k,k]) < Матем.МАШЕПС ТО
опред0:=ВКЛ;
(* пытаемся найти строку с не 0-м элементом в этом столбце *)
i:=k;
ПОКА опред0 И (i < посл) ВЫП
УВЕЛИЧИТЬ(i);
ЕСЛИ МОДУЛЬ(A[i,k]) > Матем.МАШЕПС ТО
(* переставляем эти две строки *)
ОбменСтрок(A,i,k);
опред0:=ОТКЛ;
(* соответственно меняем знак определителя *)
d:=-d
КОН
КОН
КОН;
ЕСЛИ НЕ опред0 ТО
ОТ i:=k+1 ДО посл ВЫП
ЕСЛИ МОДУЛЬ(A[i,k]) > Матем.МАШЕПС ТО
(* обнуляем k элемент этого ряда *)
t:=-A[i,k]/A[k,k];
ОТ j:=1 ДО посл ВЫП
A[i,j]:=A[i,j]+t*A[k,j]
КОН
КОН
КОН
КОН;
d:=d*A[k,k];
УВЕЛИЧИТЬ(k)
КОН;
ЕСЛИ опред0 ТО
ВОЗВРАТ 0
КОН;
ВОЗВРАТ d*A[посл,посл]
КОН Определитель;
(******************************************************************************)
ЗАДАЧА РазложитьНаLL-(A-,L+:Вид):ЦЕЛ;
(* Цель: Разложение Холецкого (A=L*L') симметрической положительно
* определенной матрицы A на две треугольные матрицы L и L', где L - это
* нижняя треугольная матрица, а L' - матрица, транспонированная к L.
******************************************************************************
* Прим: Эту задачу можно использовать для испытания на положительно
* определенность.
* До: < A > - исходная матрица A
* После: < L > - матрица L
* Ответ: 0 или НЕ_П_О *)
ПЕР
i,j,k,посл:ЦЕЛ;
сумма:Вещ;
УКАЗ
посл:=РАЗМЕР(A)-1;
ОТ k:=0 ДО посл ВЫП
сумма:=A[k,k];
ОТ j:=0 ДО k-1 ВЫП
сумма:=сумма-Матем.кв(L[k,j])
КОН;
ЕСЛИ сумма <= 0 ТО
ВОЗВРАТ НЕ_П_О
КОН;
L[k,k]:=Матем.квкор(сумма);
ОТ i:=k+1 ДО посл ВЫП
сумма:=A[i,k];
ОТ j:=0 ДО k-1 ВЫП
сумма:=сумма-L[i,j]*L[k,j]
КОН;
L[i,k]:=сумма/L[k,k]
КОН
КОН;
ВОЗВРАТ 0
КОН РазложитьНаLL;
(******************************************************************************)
ЗАДАЧА РазложитьНаLU-(A+:Вид; строки+:Перестановки):ЦЕЛ;
(* Цель: LU-разложение (A=L*U) квадратной матрицы A на два сомножителя L и U,
* где L - это нижняя треугольная матрица с единичной главной диагональю,
* а U - верхняя треугольная матрица. Эта задача используется
* совместно с РешитьИзLU для решения систем уравнений.
******************************************************************************
* До: < A > - исходная матрица
* После: < A > - содержит элементы и L и U
* <строки> - вновь созданные записи о перестановках строк
* Ответ: 0 или ВЫРОЖДЕНА *)
ПОСТ
МАЛОЕ = 1.D-20;
ПЕР
i,maxI,j,k,посл:ЦЕЛ;
вед,t,сумма:Вещ;
v:ДОСТУП К Вект.Вид;
УКАЗ
посл:=РАЗМЕР(A)-1;
СОЗДАТЬ(v,посл+1);
СОЗДАТЬ(строки,посл+1);
ОТ i:=0 ДО посл ВЫП
вед:=0;
ОТ j:=0 ДО посл ВЫП
ЕСЛИ МОДУЛЬ(A[i,j]) > вед ТО
вед:=МОДУЛЬ(A[i,j])
КОН
КОН;
ЕСЛИ вед < Матем.МАШЕПС ТО
ВОЗВРАТ ВЫРОЖДЕНА
КОН;
v[i]:=1/вед
КОН;
ОТ j:=0 ДО посл ВЫП
ОТ i:=0 ДО j-1 ВЫП
сумма:=A[i,j];
ОТ k:=0 ДО i-1 ВЫП
сумма:=сумма-A[i,k]*A[k,j]
КОН;
A[i,j]:=сумма
КОН;
maxI:=0;
вед:=0;
ОТ i:=j ДО посл ВЫП
сумма:=A[i,j];
ОТ k:=0 ДО j-1 ВЫП
сумма:=сумма-A[i,k]*A[k,j]
КОН;
A[i,j]:=сумма;
t:=v[i]*МОДУЛЬ(сумма);
ЕСЛИ t > вед ТО
вед:=t;
maxI:=i
КОН
КОН;
ЕСЛИ j # maxI ТО
ОбменСтрок(A,maxI,j);
v[maxI]:=v[j]
КОН;
строки[j]:=maxI;
ЕСЛИ A[j,j] = 0 ТО
A[j,j]:=МАЛОЕ
КОН;
ЕСЛИ j # посл ТО
t:=1/A[j,j];
ОТ i:=j+1 ДО посл ВЫП
A[i,j]:=A[i,j]*t
КОН
КОН
КОН;
ВОЗВРАТ 0
КОН РазложитьНаLU;
(******************************************************************************)
ЗАДАЧА РешитьИзLU-(A-:Вид; b-:Вект.Вид; строки:Перестановки; x+:Вект.Вид);
(* Цель: решение системы уравнений после преобразования исходной матрицы
* с помощью РазложитьНаLU
******************************************************************************
* До: < A > - матрица после РазложитьНаLU
* < b > - вектор свободных членов
* <строки> - записи о перестановках строк из РазложитьНаLU
* После: < x > - вектор решения *)
ПЕР
i,Ip,j,k,посл:ЦЕЛ;
сумма:Вещ;
УКАЗ
посл:=РАЗМЕР(A)-1;
k:=-1;
СписатьВектор(b,x);
ОТ i:=0 ДО посл ВЫП
Ip:=строки[i];
сумма:=x[Ip];
x[Ip]:=x[i];
ЕСЛИ k >= 0 ТО
ОТ j:=k ДО i-1 ВЫП
сумма:=сумма-A[i,j]*x[j]
КОН
ИНАЧЕ
ЕСЛИ сумма # 0 ТО
k:=i
КОН
КОН;
x[i]:=сумма
КОН;
ОТ i:=посл ДО 0 ПО -1 ВЫП
сумма:=x[i];
ЕСЛИ i < посл ТО
ОТ j:=i+1 ДО посл ВЫП
сумма:=сумма-A[i,j]*x[j]
КОН
КОН;
x[i]:=сумма/A[i,i]
КОН
КОН РешитьИзLU;
(******************************************************************************)
ЗАДАЧА РазложитьНаSV-(A+:Вид; s+:Вект.Вид; V+:Вид):ЦЕЛ;
(* Цель: SV разложение прямоугольной матрицы A (N x M, с N >= M) на три
* сомножителя (A=U*S*V'), где U (N x M) состоит из столбцов - собственных
* векторов, S (M x M) - диагональная матрица с элементами >= 0 (собственные
* значения) и V (M x M) ортогональная матрица. Эта задача используется
* совместно с РешитьИзSV для решения систем уравнений;
******************************************************************************
* Прим: текст этой задачи переписан из подпрограммы SVD пакета EISPACK
* До: < A > - исходная матрица
* После: < A > - содержит элементы U
* < s > - вектор собственных значений
* < V > - ортогональная матрица
* Ответ: 0 или НЕТ_СХОДИМОСТИ *)
ПЕР
i,Its,j,i1,k,l,n,посл1,посл2:ЦЕЛ;
Anorm,c,f,g,h,Scale,t,x,y,z:Вещ;
r:ДОСТУП К Вект.Вид;
label1:КЛЮЧ;
УКАЗ
посл1:=РАЗМЕР(A,0)-1;
посл2:=РАЗМЕР(A,1)-1;
g:=0;
Scale:=0;
Anorm:=0;
СОЗДАТЬ(r,посл2+1);
(* приведение к бидиагональному виду методом Хаусхолдера *)
ОТ i:=0 ДО посл2 ВЫП
l:=i+1;
r[i]:=Scale*g;
g:=0;
t:=0;
Scale:=0;
ЕСЛИ i <= посл1 ТО
ОТ k:=i ДО посл1 ВЫП
Scale:=Scale+МОДУЛЬ(A[k,i])
КОН;
ЕСЛИ Scale # 0 ТО
ОТ k:=i ДО посл1 ВЫП
A[k,i]:=A[k,i]/Scale;
t:=t+Матем.кв(A[k,i])
КОН;
f:=A[i,i];
g:=-Матем.знак(f)*Матем.квкор(t);
h:=f*g-t;
A[i,i]:=f-g;
ЕСЛИ i # посл2 ТО
ОТ j:=l ДО посл2 ВЫП
t:=0;
ОТ k:=i ДО посл1 ВЫП
t:=t+A[k,i]*A[k,j]
КОН;
f:=t/h;
ОТ k:=i ДО посл1 ВЫП
A[k,j]:=A[k,j]+f*A[k,i]
КОН
КОН
КОН;
ОТ k:=i ДО посл1 ВЫП
A[k,i]:=Scale*A[k,i]
КОН
КОН
КОН;
s[i]:=Scale*g;
g:=0;
t:=0;
Scale:=0;
ЕСЛИ (i <= посл1) И (i # посл2) ТО
ОТ k:=l ДО посл2 ВЫП
Scale:=Scale+МОДУЛЬ(A[i,k])
КОН;
ЕСЛИ Scale # 0 ТО
ОТ k:=l ДО посл2 ВЫП
A[i,k]:=A[i,k]/Scale;
t:=t+Матем.кв(A[i,k])
КОН;
f:=A[i,l];
g:=-Матем.знак(f)*Матем.квкор(t);
h:=f*g-t;
A[i,l]:=f-g;
ОТ k:=l ДО посл2 ВЫП
r[k]:=A[i,k]/h
КОН;
ЕСЛИ i # посл1 ТО
ОТ j:=l ДО посл1 ВЫП
t:=0;
ОТ k:=l ДО посл2 ВЫП
t:=t+A[j,k]*A[i,k]
КОН;
ОТ k:=l ДО посл2 ВЫП
A[j,k]:=A[j,k]+t*r[k]
КОН;
КОН
КОН;
ОТ k:=l ДО посл2 ВЫП
A[i,k]:=Scale*A[i,k]
КОН
КОН
КОН;
Anorm:=Матем.макс(Anorm,МОДУЛЬ(s[i])+МОДУЛЬ(r[i]));
КОН;
(* накопление правых преобразований *)
ОТ i:=посл2 ДО 0 ПО -1 ВЫП
ЕСЛИ i < посл2 ТО
ЕСЛИ g # 0 ТО
ОТ j:=l ДО посл2 ВЫП
(* двойное деление предохраняет от возможной потери точности *)
V[j,i]:=(A[i,j]/A[i,l])/g
КОН;
ОТ j:=l ДО посл2 ВЫП
t:=0;
ОТ k:=l ДО посл2 ВЫП
t:=t+A[i,k]*V[k,j]
КОН;
ОТ k:=l ДО посл2 ВЫП
V[k,j]:=V[k,j]+t*V[k,i]
КОН
КОН
КОН;
ОТ j:=l ДО посл2 ВЫП
V[i,j]:=0;
V[j,i]:=0
КОН
КОН;
V[i,i]:=1;
g:=r[i];
l:=i
КОН;
(* накопление левых преобразований *)
ОТ i:=посл2 ДО 0 ПО -1 ВЫП
l:=i+1;
g:=s[i];
ЕСЛИ i < посл2 ТО
ОТ j:=l ДО посл2 ВЫП
A[i,j]:=0
КОН
КОН;
ЕСЛИ g # 0 ТО
ЕСЛИ i # посл2 ТО
ОТ j:=l ДО посл2 ВЫП
t:=0;
ОТ k:=l ДО посл1 ВЫП
t:=t+A[k,i]*A[k,j]
КОН;
(* двойное деление предохраняет от возможной потери точности *)
f:=(t/A[i,i])/g;
ОТ k:=i ДО посл1 ВЫП
A[k,j]:=A[k,j]+f*A[k,i]
КОН;
КОН
КОН;
ОТ j:=i ДО посл1 ВЫП
A[j,i]:=A[j,i]/g
КОН;
ИНАЧЕ
ОТ j:=i ДО посл1 ВЫП
A[j,i]:=0
КОН
КОН;
A[i,i]:=A[i,i]+1
КОН;
(* диагонализация бидиагональной формы *)
ОТ k:=посл2 ДО 0 ПО -1 ВЫП
Its:=0;
КОЛЬЦО
l:=k;
(* проверка на расщепление *)
КОЛЬЦО
ЕСЛИ l < 0 ТО
label1:=ВКЛ;
ВЫХОД
КОН;
n:=l-1;
ЕСЛИ (МОДУЛЬ(r[l])+Anorm) = Anorm ТО
label1:=ОТКЛ;
ВЫХОД
КОН;
(* r[0] всегда =0, поэтому s[-1] никогда не будет вычисляться *)
ЕСЛИ (МОДУЛЬ(s[n])+Anorm) = Anorm ТО
label1:=ВКЛ;
ВЫХОД
КОН;
УМЕНЬШИТЬ(l)
КОН; (* кольцо *)
ЕСЛИ label1 ТО
c:=0;
t:=1;
i:=l;
КОЛЬЦО
ЕСЛИ i > k ТО
ВЫХОД
КОН;
f:=t*r[i];
r[i]:=c*r[i];
ЕСЛИ (МОДУЛЬ(f)+Anorm) = Anorm ТО
ВЫХОД
КОН;
g:=s[i];
h:=Матем.Пифагор(f,g);
s[i]:=h;
c:=g/h;
t:=-f/h;
ОТ j:=0 ДО посл1 ВЫП
y:=A[j,n];
z:=A[j,i];
A[j,n]:=y*c+z*t;
A[j,i]:=-y*t+z*c
КОН;
УВЕЛИЧИТЬ(i)
КОН (* кольцо *)
КОН;
(* проверка сходимости *)
z:=s[k];
ЕСЛИ l = k ТО
ЕСЛИ z < 0 ТО
(* сделаем s[k] неотрицательным *)
s[k]:=-z;
ОТ j:=0 ДО посл2 ВЫП
V[j,k]:=-V[j,k]
КОН
КОН;
(* сходимость *)
ВЫХОД
КОН;
ЕСЛИ Its = 30 ТО
ВОЗВРАТ НЕТ_СХОДИМОСТИ
КОН;
УВЕЛИЧИТЬ(Its);
(* сдвиг снизу на минор 2 х 2 *)
x:=s[l];
n:=k-1;
y:=s[n];
g:=r[n];
h:=r[k];
f:=0.5*(((g+z)/h)*((g-z)/y)+y/h-h/y);
g:=Матем.Пифагор(f,1);
f:=x-(z/x)*z+(h/x)*(y/(f+Матем.знак(f)*МОДУЛЬ(g))-h);
(* следующее QR преобразование *)
c:=1;
t:=1;
ОТ i1:=l ДО n ВЫП
i:=i1+1;
g:=r[i];
y:=s[i];
h:=t*g;
g:=c*g;
z:=Матем.Пифагор(f,h);
r[i1]:=z;
c:=f/z;
t:=h/z;
f:=x*c+g*t;
g:=-x*t+g*c;
h:=y*t;
y:=y*c;
ОТ j:=0 ДО посл2 ВЫП
x:=V[j,i1];
z:=V[j,i];
V[j,i1]:=x*c+z*t;
V[j,i]:=-x*t+z*c
КОН;
z:=Матем.Пифагор(f,h);
s[i1]:=z;
(* вращение может быть произвольным при z=0 *)
ЕСЛИ z # 0 ТО
c:=f/z;
t:=h/z
КОН;
f:=c*g+t*y;
x:=-t*g+c*y;
ОТ j:=0 ДО посл1 ВЫП
y:=A[j,i1];
z:=A[j,i];
A[j,i1]:=y*c+z*t;
A[j,i]:=-y*t+z*c
КОН
КОН;
r[l]:=0;
r[k]:=f;
s[k]:=x
КОН (* кольцо *)
КОН;
ВОЗВРАТ 0
КОН РазложитьНаSV;
(******************************************************************************)
ЗАДАЧА ОбнулитьS-(s+:Вект.Вид; откл:Вещ);
(* Цель: обнуляет наименьшие собственные значения по заданной границе
******************************************************************************
* До: < s > - вектор собственных значений S
* <откл> - относительное отклонение (граница будет равна откл*макс(S) )
* После: < s > - измененный вектор собственных значений *)
ПЕР
граница:Вещ;
i:ЦЕЛ;
УКАЗ
граница:=откл*Вект.макс(s);
ОТ i:=0 ДО РАЗМЕР(s)-1 ВЫП
ЕСЛИ s[i] < граница ТО
s[i]:=0
КОН
КОН
КОН ОбнулитьS;
(******************************************************************************)
ЗАДАЧА РешитьИзSV-(U-:Вид; s-:Вект.Вид; V-:Вид; b-,x+:Вект.Вид);
(* Цель: решение системы уравнений после преобразования исходной матрицы
* с помощью РазложитьНаSV и обнуления наименьших собственных значений
* с помощью ОбнулитьS
******************************************************************************
* До: < U > - матрица U после РазложитьНаSV
* < s > - вектор S после РазложитьНаSV
* < V > - матрица V после РазложитьНаSV
* < b > - вектор свободных членов
* После: < x > - вектор решения (V*Diag(1/S(i))*U'*b, для S(i) # 0) *)
ПЕР
i,j,посл1,посл2:ЦЕЛ;
сумма:Вещ;
врем:ДОСТУП К Вект.Вид;
УКАЗ
посл1:=РАЗМЕР(U,0)-1;
посл2:=РАЗМЕР(U,1)-1;
СОЗДАТЬ(врем,посл2+1);
ОТ j:=0 ДО посл2 ВЫП
сумма:=0;
ЕСЛИ s[j] > 0 ТО
ОТ i:=0 ДО посл1 ВЫП
сумма:=сумма+U[i,j]*b[i]
КОН;
сумма:=сумма/s[j]
КОН;
врем[j]:=сумма
КОН;
ОТ j:=0 ДО посл2 ВЫП
сумма:=0;
ОТ i:=0 ДО посл2 ВЫП
сумма:=сумма+V[j,i]*врем[i]
КОН;
x[j]:=сумма
КОН;
КОН РешитьИзSV;
(******************************************************************************)
ЗАДАЧА ВосстановитьИзSV-(U-:Вид; s-:Вект.Вид; V-,A+:Вид);
(* Цель: восстанавливает матрицу A перемножая U, S и V' после
* обнуления наименьших собственных значений задачей ОбнулитьS
******************************************************************************
* До: < U > - матрица U после РазложитьНаSV
* < s > - вектор S после РазложитьНаSV
* < V > - матрица V после РазложитьНаSV
* После: < A > - восстановленная матрица *)
ПЕР
i,j,k,посл1,посл2:ЦЕЛ;
УКАЗ
посл1:=РАЗМЕР(U,0)-1;
посл2:=РАЗМЕР(U,1)-1;
ОТ i:=0 ДО посл1 ВЫП
ОТ j:=0 ДО посл2 ВЫП
A[i,j]:=0;
ОТ k:=0 ДО посл2 ВЫП
ЕСЛИ s[k] > 0 ТО
A[i,j]:=A[i,j]+U[i,k]*V[j,k]
КОН
КОН
КОН
КОН
КОН ВосстановитьИзSV;
(******************************************************************************)
ЗАДАЧА РазложитьНаQR-(A+,R+:Вид):ЦЕЛ;
(* Цель: QR разложение прямоугольной матрицы (N x M, с N >= M) на два
* сомножителя (A=Q*R), где Q (N x M) ортогональная матрица, а
* R (M x M) - верхняя треугольная матрица. Эта задача используется
* совместно с РешитьИзQR для решения систем уравнений.
******************************************************************************
* До: < A > - исходная матрица
* После: < A > - содержит элементы Q
* < R > - верхняя треугольная матрица
* Ответ: 0 или ВЫРОЖДЕНА *)
ПЕР
i,j,k,посл1,посл2:ЦЕЛ;
сумма:Вещ;
УКАЗ
посл1:=РАЗМЕР(A,0)-1;
посл2:=РАЗМЕР(A,1)-1;
ОТ k:=0 ДО посл2 ВЫП
(* вычисляем k-й диагональный элемент *)
сумма:=0;
ОТ i:=0 ДО посл1 ВЫП
сумма:=сумма+Матем.кв(A[i,k])
КОН;
ЕСЛИ сумма = 0 ТО
ВОЗВРАТ ВЫРОЖДЕНА
КОН;
R[k,k]:=Матем.квкор(сумма);
(* делим k-й столбец < A > на k-й диагональный элемент
* (это начало перезаписи < A > на < q >) *)
ОТ i:=0 ДО посл1 ВЫП
A[i,k]:=A[i,k]/R[k,k]
КОН;
ОТ j:=(k+1) ДО посл2 ВЫП
(* дополняем оставшиеся элементы в столбце *)
сумма:=0;
ОТ i:=0 ДО посл1 ВЫП
сумма:=сумма+A[i,k]*A[i,j]
КОН;
R[k,j]:=сумма;
(* обновляем элементы столбца матрицы < q > (< A >) *)
ОТ i:=0 ДО посл1 ВЫП
A[i,j]:=A[i,j]-A[i,k]*R[k,j]
КОН
КОН
КОН;
ВОЗВРАТ 0;
КОН РазложитьНаQR;
(******************************************************************************)
ЗАДАЧА РешитьИзQR-(q-,R-:Вид; b-,x+:Вект.Вид);
(* Цель: решение системы уравнений после преобразования исходной матрицы
* с помощью РазложитьНаQR
******************************************************************************
* До: < q > - матрица Q после РазложитьНаQR
* < R > - матрица R после РазложитьНаQR
* < b > - вектор свободных членов
* После: < x > - вектор решения *)
ПЕР
i,j,посл1,посл2:ЦЕЛ;
сумма:Вещ;
УКАЗ
посл1:=РАЗМЕР(q,0)-1;
посл2:=РАЗМЕР(q,1)-1;
(* строим q'b и записываем результат в x *)
ОТ j:=0 ДО посл2 ВЫП
x[j]:=0;
ОТ i:=0 ДО посл1 ВЫП
x[j]:=x[j]+q[i,j]*b[i]
КОН
КОН;
(* обновляем вектором решения *)
x[посл2]:=x[посл2]/R[посл2,посл2];
ОТ i:=(посл2-1) ДО 0 ПО -1 ВЫП
сумма:=0;
ОТ j:=(i+1) ДО посл2 ВЫП
сумма:=сумма+R[i,j]*x[j]
КОН;
x[i]:=(x[i]-сумма)/R[i,i]
КОН
КОН РешитьИзQR;
КОН Матр.
|
|
Программируем по-русски